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Summary 


Background 

On the 17 th of January 2014, the minister of Economie Affairs decided to reduce 
production from five production clusters in the center of the Groningen field to 3 
bem per year for the period 2014 - 2016 in order to try to reduce the seismicity in 
the center of the field. In addition, total field production was limited to 42.5 bem for 
2014 and 2015 and 40 bem for 2016 (EZ 2014). Preceding this decision, technical 
reports (NAM 2013, TNO 2013) concluded that the seismicity is related to the 
compaction (and hence the production) of the Groningen field. 

In 2014 42,41 bem was produced from the Groningen gas field. Additional reports 
(TNO 2014a, TNO 2014b) indicated a change in the event density of the field. A 
decrease in the center of the field was shown as well as slight increases north of 
Hoogezand and nearby Tjuchem. 

Scope 

In the beginning of 2015, the Minister of Economie Affairs decided to impose 
production caps to the Groningen gas production on a semi-yearly basis. For the 
first six months of 2015 a maximum gas production of 16,5 bem was allowed. In 
January 2015 State Supervision of Mines (SSM) advised to reduce production to 
39,4 bem in 2015 and 2016. After a new advice from SSM, which is expected to be 
presented in June 2015, the definitive production maximum for 2015 will be defined. 

In support of their advice of June 1 st 2015, State Supervision of Mines has 
requested the following additional technical evaluations from TNO-AGE: 

• An update on the seismicity of the Groningen field 

• compaction field based on inversion of subsidence data 

Gas production in 2014 

In 2014 the total gas production of the Groningen gas field was 42,41 bem, which is 
less than the imposed production cap of 42,5 bem. The Loppersum production 
clusters (LRM, PAU, POS, OVS, and ZND) produced in total 2,57 bem in 2014, 
which is below the production cap of 3 bem. Production varied over the year, with 
the majority of gas being produced during the winter months. 

Update on the seismicity of the Groningen field 2014/2015 

The distribution of higher magnitude (M L >2) events occurring since September 2014 

can be explained by 

The distances to producing clusters vs non-producing clusters which can 
explain the events close to Appingedam 

The increase of production at the Ten Post cluster (POS) in December 
2014 

The first explanation would indicate that the effect of reducing production of the 
Loppersum clusters has been partially overruled by production of other clusters 
close to Appingedam. The second explanation would indicate that sudden 
increases in production could lead to a changing pattern of events in time and 
space. The latter statement cannot yet be proven with statistical significance and 
should therefore be further evaluated. An analysis of production and seismic events 
occurring over time could possibly provide further statistical significance. 
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The observed density of seismic events from April 2014 to April 2015 is different 
from densities observed during earlier years. Largest seismic event densities are 
concentrated in the Southwest while the center of the field is marked by lower 
densities. In previous years the density of seismic events in the center were highest 
in the Groningen field (see also TNO 2014b). This indicates that the reduction of 
production in the central area has a marked influence on the number of events in 
the same area. Additionally, there is a striking match of the event density in 
2014/2015 to two known fault systems in the field. These fault systems correspond 
to areas in the field where differential compaction, known to be an indicator tor the 
occurrence of seismic events, exists (Figure ii). 

Statistical analysis on the number of seismic events indicates that the number of 
events per day in the center of the field has halved since January 2014. The 
southwestern area, however, shows an increase in the number of events per day. 
Similar to the production, the seismic events of the Groningen fields exhibits clear 
seasonality with a lag of some two months between production changes and a 
change in seismic events. 

A Bayesian change point model that has been successfully applied in Oklahoma, 
U.S.A., has also been applied to the Groningen gas field. A change point in the 
center of the field is found in January 2003. Event rates after 2003 have quadrupled 
compared to the years prior to 2003. This would indicate that the fault system in the 
center of the Groningen field has reached criticality in the beginning of 2003, i.e. 
small changes in stress will lead to seismic events. Over the whole of the field 
change points are identified which vary in time (from 2003 to 2010); the earliest 
times in the center of the field and later times at the edges of the field. This 
corresponds to the observation that events have started to occur in the center of the 
field and have spread in time over the field. If the change point indicates when a 
fault system becomes critical then this also means that different fault systems have 
become critical at different times. 

After 2009 no change point in the center of the field has been found for seismic 
events with magnitudes larger than M L =1.5. The number of events (ML>1.5) since 
2014 is probably not enough to show a change point in event rates. Thus this data 
cannot be used to prove statistically significant changes in event rates since the 
production reduction of January 2014. 
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Figure i. Event density (number of events per km 2 ) from April 1 st 2014 to April 1 st 2015 shown with 
the faults in the reservoir (dark red) and the contour of the field (dark blue). 

Alternative compaction field 

Inversion of subsidence data has provided a correction to the compaction field 
presented in TNO (2013, 2014a,b). The correction is predominantly applied in 
regions where previously erroneous porosity estimations or aquifer activity were 
suspected. The area of maximum compaction has shifted to the west and does not 
correspond to the area of maximum event density in the center of the field. This 
indicates that in this regard the presence of faults is more important tor seismicity 
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than the compaction itself. Also differential compaction, known to be an indicator tor 
the occurrence of seismicity, is visible over faults. 

This leads to the conclusion that the existing seismological model which NAM has 
used in the production plan (NAM, 2013), based on an empirical relation between 
total compaction and the occurrence of events, needs to be updated. As indicated 
in TNO (2013, 2014a,b) the faults in the reservoir play an important role in the 
occurrence of events within the field and therefore they have to be taken into 
account in any future seismological model. 


228000 232000 236000 240000 244000 248000 2S2000 256000 260000 264000 268000 
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Figure ii. Compaction (m) in 2013 obtained through inversion of subsidence measurements. The 
red line shows the contour of the Groningen field and the black lines are the faults that are present 
in the geological model in Petrel (NAM, 2013). Also shown is the seismicity in the field, the size of 
the symbols indicates the magnitudes of the events 
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Figures 

Figure 2-1. OverView of the mentioned clusters: the Loppersum clusters (LRM, 

PAU, POS, OVS, ZND); the clusters close to Appingedam (AMR, SDB, TJM) 
and Eemskanaal (EKL). Additionally a few cities are indicated 
(GRO=Groningen, HGZ= Hoogezand, WIN=Winschoten, DLZ= Delfzijl, 
LOP=Loppersum, APD=Appingedam). The contour of the Groningen gas field 
and the Annerveen gas field to the south is indicated in green, the topography 

in black.11 

Figure 2-2. Production (taken from www.nlog.nl) of the five Loppersum clusters: 
Leermens (LRM), Overschild (OVS), De Paauwen (PAU), Ten Post (POS), ‘t 

Zandt (ZND) in 2014 indicated per month.12 

Figure 2-3. Event density (number of events per km 2 ) from April 1 st 2014 to April 1 st 
2015. The observed events and their magnitudes are indicated by the colored 

small circles.15 

Figure 2-4. Event density (number of events per km 2 ) from April 1 st 2014 to April 1 st 
2015 shown with the faults in the reservoir (dark red) and the contour of the 

field (dark blue).16 

Figure 2-5. Event density (number of events per km 2 ) from April 1 st 2012 to April 1 st 

2013. The observed events and their magnitudes are indicated by the colored 

small circles.17 

Figure 2-6. Event density (number of events per km 2 ) from April 1 st 2013 to April 1 st 

2014. The observed events and their magnitudes are indicated by the colored 

small circles.18 

Figure 2-7. Difference in event density (number of events per km 2 ) between April 1 st 
2014 - April 1 st 2015 (Figure 2-3) and April 1 st 2013-April 1 st 2014 (Figure 2-6) 

A negative (green) difference indicates a lower event density in 2014/2015 

compared to 2013/2014.19 

Figure 2-8. Difference in event density (number of events per km 2 ) between April 1 st 
2014 - April 1 st 2015 (Figure 2-3) and April 1 st 2012-April 1 st 2013 (Figure 2-5) 

A negative (green) difference indicates a lower event density in 2014/2015 

compared to 2012/2013.20 

Figure 2-9. Number of events occurring within the contour of the Groningen gas 

field as a function of time and Magnitude (M).23 

Figure 2-10. Autocorrelation of the production on a monthly basis.23 

Figure 2-11. The correlation between the production on a monthly basis and the 

number of seismic events.24 

Figure 2-12. The correlation between the production on a monthly basis and the 

change in seismic events on a monthly basis.24 

Figure 2-13. Autocorrelation of the production on a monthly basis.25 

Figure 2-14a.) The correlation between the production on a monthly basis and the 
number of seismic events and b.) The correlation between the production on a 

monthly basis and the change in seismic events on a monthly basis.25 

Figure 2-15. The probability of change in time over the period of 1991 up to now. ..27 
Figure 2-16a.) The pre change date event rate (in events/day) and b.) the post 
change date event rate (in events/day) 
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Figure 2-17. Time of event rate changes evaluated at 50 local points in the 

Groningen field.28 

Figure 3-1. Schematics of the forward method.29 

Figure 3-2 Prior compaction fields (top row) and estimated compaction fields 

(bottom row) in 1993 (left) and in 2013 (right).30 

Figure 3-3. Compaction (m) in 2013 obtained through inversion of subsidence 
measurements (section 3.1). The red line gives the contour of the Groningen 
field and the black lines are the faults that are present in the geological model 

in Petrel (NAM, 2013).31 

Figure 3-4. Compaction (m) in 2013 obtained through inversion of subsidence 


measurements (section 3.1). The red line gives the contour of the Groningen 
field and the black lines are the faults that are present in the geological model 
in Petrel (NAM, 2013). Also shown is the seismicity in the field, the size of the 

symbols indicates the magnitudes of the events.32 

Figure 3-5. Compaction (m) in 2013 obtained through inversion of subsidence 
measurements (section 3.1). The red line gives the contour of the Groningen 


field and the faults in the geological Petrel model (NAM, 2013) are indicated 

with their offset.33 

Figure 3-6. The difference between the compaction field of TNO (2013, 2014a,b) 

and the compaction field resulting from inversion (m).34 

Figure 3-7. Figure 5.13 from TNO (2013). Compaction in 2012 calculated with the 
RTiCM model using the subsurface model.35 
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Tables 

Table 2-1. Production in 2014 of the five Loppersum clusters.12 

Table 2-2. Induced seismicity (taken trom www.knmi.nl) of the Groningen field, 

events larger than M L =2 and after September 2014.13 

Table 2-3. The number of events in the regions Central, SW and Other as a function 
of the number of days since the start of observed seismicity on December 5 th 

1991.21 

Table 2-4. The event rate, including Standard deviation, in the regions Central, SW 
and Other as a function of the number of days since the start of seismicity on 
December 5 th 1991.22 
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1 Introduction 


Background 

On the 17 th of January 2014, the minister of Economie Affairs decided to reduce 
production from five production clusters in the center of the Groningen field to 3 
bem per year for the period 2014 - 2016 in order to try to reduce the seismicity in 
the center of the field. In addition, total field production was limited to 42.5 bem for 
2014 and 2015 and 40 bem for 2016 (EZ 2014). Preceding this decision, technical 
reports (NAM 2013, TNO 2013) concluded that the seismicity is related to the 
compaction (and hence the production) of the Groningen field. 

In 2014 42,41 bem was produced from the Groningen gas field. Additional reports 
(TNO 2014a, TNO 2014b) indicated a change in the rate of seismicity in the field. A 
decrease in the center of the field was shown as well as slight increases north of 
Hoogezand (Southwest area) and nearby Tjuchem (eastern area). 

Scope 

In the beginning of 2015, the Minister of Economie Affairs decided to impose 
production caps to the Groningen gas production on a semi-yearly basis. For the 
first six months of 2015 a maximum gas production of 16,5 bem was allowed. In 
January 2015 State Supervision of Mines (SSM) advised to reduce production to 
39,4 bem in 2015 and 2016. After a new advice from SSM, which is expected to be 
presented in June 2015, the definitive production maximum for 2015 will be defined. 

In support of their advice, State Supervision of Mines has requested the following 
additional technical evaluations from TNO-AGE: 

• An update on the seismicity of the Groningen field 

• compaction field based on inversion of subsidence data 


Report setup 

Chapters 2 and 3 report the key results and findings of TNO’s evaluations. 

In chapter 2 the seismicity of the Groningen field since January 2014 is reported 
including statistical analysis on observed seismic events. Chapter 3 presents the 
compaction field I based on the inversion of subsidence data and its implications for 
the link between compaction, seismicity and the existing faults in the reservoir. 
Finally, chapter 4 summarizes the findings from chapters 2 and 3 with regards to 
the questions of SSM. 
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2 Recent developments of the events in 2014 and 
2015 


2.1 Introduction 

On January 17 ,h 2014 a reduction of production was imposed on the Groningen gas 
field. As a part of this overall reduction, the production in the five clusters in the 
center of the field (clusters Leermens (LRM), Ten Post (POS), de Paauwen (PAU), 
Overschild (OVS), ‘t Zandt (ZND); see Figure 2-1) was reduced to 3 bcm per year. 

TNO (2013) and NAM (2013) concluded that seismicity in the field is linked to 
compaction, and compaction on its turn is directly linked to production. Hence the 
observed rate of seismic events (number of events per unit time and unit area) 
should be indicative for whether the decrease in production has had an effect on 
seismicity since the implementation of the production reduction measures. 



x [m] x 10 5 


Figure 2-1. OverView of the mentioned clusters: the Loppersum clusters (LRM, PAU, POS, OVS, 
ZND); the clusters close to Appingedam (AMR, SDB, TJM) and Eemskanaal (EKL). 
Additionally a few cities are indicated (GRO=Groningen, HGZ= Hoogezand, 
WIN=Winschoten, DLZ= Delfzijl, LOP=Loppersum, APD=Appingedam). The contour of 
the Groningen gas field and the Annerveen gas field to the south is indicated in green, 
the topography in black. 
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This chapter presents an update on chapter 2 of TNO (2014b) regarding the event 
densities. In section 2.2 the production of the field in 2014 is presented, in section 
2.3 the observed seismic events in 2014 and 2015 above magnitude M L =2.0 are 
discussed. In section 2.4 the event density from April 1 st 2014 to April 1 st 2015 is 
shown. The event densities are compared to the observed event densities for the 
same period in the previous years in section 2.5. Section 2.6 describes the statistics 
of the induced seismicity in a Bayesian analysis and an analysis for seasonality of 
the events. Finally in section 2.7 a Bayesian Point Change model is applied to the 
induced events. 


2.2 Gas Production of the Groningen field 

With 42,41 bcm of gas produced in 2014, the total production of the Groningen field 
stayed below the imposed production cap of 42,5 bcm. The so-called Loppersum 
clusters (LRM, PAU, POS, OVS, ZND) produced in total 2,57 bcm in 2014 (see 
Table 2-1), which is also below the imposed production cap of 3 bcm/yr. Production 
varied over the year and was highest during the winter months (Figure 2-2). 

Table 2-1. Production in 2014 of the five Loppersum clusters. 


Cluster 

Production in 2014 
(bcm) 

LRM 

0.57 

PAU 

0.28 

POS 

0.61 

OVS 

0.62 

ZND 

0.49 



Figure 2-2. Production (taken from www.nlog.nl) of the five Loppersum clusters: Leermens (LRM), 
Overschild (OVS), De Paauwen (PAU), Ten Post (POS), ‘t Zandt (ZND) in 2014 
indicated per month. 
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2.3 Induced Seismicity of the Groningen field in 2014/2015, events larger than 
M l =2 


Table 2-2 shows the events with magnitudes larger than M L =2 that have occurred in 
Groningen since September 2014. This date has been chosen as the production 
reduction in Loppersum will have influenced part of the field in the center (TNO 
2014a,b), influencing the number and possibly magnitude of the events. The 
underlying assumption here is the pressure diffusion model, which was presented in 
TNO (2014b). 

Table 2-2. Induced seismicity (taken trom www.knmi.nl) of the Groningen field, events larger than 
Ml=2 and after September 2014. 


Event 

date 


Garnerwolde 

30-09-2014 

2.8 

Zandeweer 

05-11-2014 

2.9 

Woudbloem 

30-12-2014 

2.8 

Wirdum 

06-01-2015 

2.7 

Appingedam 

25-02-2015 

2.3 

Appingedam 

24-03-2015 

2.3 


The events in Garnerwolde and Woudbloem occurred in the Southwest region of the 
field, not affected by the production reduction in the center. In TNO (2014b) the 
Southwest region is described in more detail. 

The Zandeweer event occurred in the north of the field. This part of the field has not 
yet been influenced by the production reduction of the Loppersum clusters as the 
travel speed of the pressure is influenced by the permeability of the reservoir (for 
details see TNO, 2014b). 

The two events near Appinqedam both occurred in 2015. These events occurred 
close to the production clusters Amsweer (AMR), Siddeburen (SDB) and Tjuchem 
(TJM) as well as to the Loppersum clusters Overschild (OVS) and Leermens (LRM) 
(Figure 2-1). The increase of events near Appingedam may indicate that the 
pressure wave associated with the continuing production in the nearby clusters of 
AMR, SDB and TJM causes compaction and consequently seismicity in that area. It 
is, however, too early to draw conclusions from this statement and more 
observations are needed to support statistical significance. 

The event near Wirdum occurred in the Loppersum area in the beginning of 
January 2015. With a distance of 2.1 km the POS cluster is the nearest cluster. The 
sudden increase of production at this cluster in December 2014 (Figure 2-2), may 
have induced this event. In this case the pressure wave would have traveled 
between 1.4 km and 2.5 km in one month, depending on the reservoir permeability 
(150 - 500 mD). Again, more observations from induced events are required to 
support statistical significance to substantiate such a direct relation between 
seismic events and a sudden increase in production. 
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2.4 Observed event density from April 2014 to April 2015 

Figure 2-3 shows the observed event density tor the period from April 1 st 2014 to 
April 1 st 2015. The event density was determined using a Kernei Density (Standard 
GIS application) with a radius of 5 km and a cell size of 50 m. As indicated in TNO 
(2014b), the pressure wave should have traveled approximately 2 to 4 km between 
January 17 th 2014 and April 1 st 2014. Therefore April 1 st 2014 is chosen as the date 
after which the rate of seismic events will possibly be affected by the reduction of 
production. 

The average event density is around 0,25 events per km 2 with largest densities in 
the Southwest periphery of the field. Other areas with increased (with respect to the 
background) event densities during this period are 1) the area west of Delfzijl 
(Appingedam), 2) the area to the north between Middelstum and Loppersum and 3) 
the area near Tjuchem (Figure 2-3). Compared to TNO (2014b; Figure 2-5) the 
areas marked by high event densities correspond well except for the Appingedam 
area which appears to be characterized by a higher event density now, as shown in 
Figure 2-3. This means that in the period from November 1 st 2014 to April 1 st 2015 
events have occurred in the Appingedam area. As mentioned previously, the effects 
of the reduction of production at the nearby Leermens (LRM) cluster to the 
northwest of Appingedam and the Overschild (OVS) cluster to the Southwest of 
Appingedam seem to be overruled by the ongoing production to the southeast and 
south of Appingedam (Figure 2-1). 

Figure 2-4 shows the event density for the period of April 1 st 2014 to April 1 st 2015 
together with the faults in the reservoir. The match of the event density to two 
known, mainly NW-SE trending faults Systems in the field is striking. One active 
fault system in the north of the field stretches from the northwest to the east and 
another active fault system is located in the Southwest of the field. 
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Figure 2-3. Event density (number of events per km 2 ) from April 1 st 2014 to April 1 st 2015. The 
observed events and their magnitudes are indicated by the colored small circles. 
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Figure 2-4. Event density (number of events per km 2 ) from April 1 st 2014 to April 1 st 2015 shown 
with the faults in the reservoir (dark red) and the contour of the field (dark blue). 

2.5 Comparison to earlier years (2012 and 2013) 

In Figure 2-5 and Figure 2-6 the event density is shown tor the period April 1 st 2012 
to April 1 st 2013 (Figure 2-5) and April 1 st 2013 to April 1 st 2014 (Figure 2-6). 
Compared to Figure 2-3 the amplitudes of the event densities are larger by a factor 
of 2 (up to 0.5 events per km 2 ). Also the shape of the event density is different. The 
largest event densities are observed in the center of the field and in the 
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Eemskanaal area (Figure 2-1). In Figure 2-7 and Figure 2-8 the difference between 
the event density in 2014/2015 and the event density in 2013/2014 (Figure 2-7) and 
the event density in 2012/2013 (Figure 2-8) is shown. In the center a clear decrease 
in the event density is visible tor both difference maps. Thus the event density in the 
center of the field has diminished since the production reduction of January 2014 
compared to previous years. 



Earthquakes between April Ist 2012 to April Ist 2013 
® <-2 
© 2-3 

(O)»* 

B 3roninger gasveld 
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Figure 2-5. Event density (number of events per km 2 ) trom April 1 st 2012 to April 1 st 2013. The 
observed events and their magnitudes are indicated by the colored small circles. 
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A 



Figure 2-6. Event density (number of events per km 2 ) from April 1 st 2013 to April 1 st 2014. The 
observed events and their magnitudes are indicated by the colored small circles. 
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Figure 2-7. Difference in event density (number of events per km 2 ) between April 1 st 2014 - April 1 st 
2015 (Figure 2-3) and April 1 st 2013-April 1 st 2014 (Figure 2-6) A negative (green) 
difference indicates a lower event density in 2014/2015 compared to 2013/2014. 
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Figure 2-8. Difference in event density (number of events per km 2 ) between April 1 st 2014 - April 1 st 
2015 (Figure 2-3) and April 1 st 2012-April 1 st 2013 (Figure 2-5) A negative (green) 
difference indicates a lower event density in 2014/2015 compared to 2012/2013. 

2.6 Statistics of the induced seismicity of the Groningen field 

In TNO (2014b) a Bayesian analysis of the event rate from 1991 until November 
2014 is presented. Section 2.6.1 provides an update to this analysis. The driver for 
these statistical analyses was to assess whether a significant change in the 
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occurrence of events has developed since the production reduction in the center of 
the field. In this section two different Bayesian methods are shown. 

2.6.1 Previous Bayesian analysis 

The Bayesian analysis in the report of TNO (2014b) assumed three basic trend 
models: a Poison distribution with a constant, an increasing, and a decreasing 
seismic event rate. In the period 2003 - January 17 th 2014 the rates in three areas 
of the Groningen field were shown to be increasing. The main question was to 
investigate whether this trend of increasing events was halted. It was concluded 
that the chosen method slightly favored a decreasing model, but the data were too 
sparse to come to a firm conclusion. It turns out that this method will only give a 
statistically significant answer in 5 to 10 years. Therefore an alternative was used in 
the next section. 

2.6.2 Alternative Bayesian analysis 

In order to strengthen the results of the previous study we have chosen an 
alternative approach here. We subdivided the data in segments of 1000 days and 
looked at a constant rate model for each segment. The Poisson model for the rate a 
is then given as 

p (k | a) = (aT) k exp (-aT) / k! 

In each segment we determined the mode of a of the posterior distribution p(a | k), 
assuming a constant prior p(a). The maximum value of the posterior distribution in 
each segment is reached for the value a = k / T. 


Table 2-3. The number of events in the regions Central, SW and Other as a function of the number 
of days since the start of observed seismicity on December 5 th 1991. 


Time (days) 

Central 

Events M L £1 

SW 

Other 

0-1000 

7 

2 

20 

1000-2000 

7 

0 

14 

2000 - 3000 

12 

5 

17 

3000 - 4000 

7 

3 

11 

4000 - 5000 

29 

7 

23 

5000 - 6000 

31 

12 

28 

6000 - 7000 

31 

10 

49 

7000 - 8080 

63 

19 

106 

8080 - 8539 

11 

18 

42 


The database of the KNMI is used to evaluate the number of events in a given time 
period and region of the Groningen field. All magnitudes above M L =1 have been 
taken into account since the first event on December 5 th 1991. The magnitude of 
completeness (which is the magnitude from which all events over the Groningen 
field have been registered) changes over the period 1991-2015. The seismometer 
network was significantly extended up to 1996. The first two time periods (0 - 2000 
days) will therefore not have all seismic events included in the database. The 
magnitude of completeness from 1996 over the whole of the Groningen field is 
M l =1 .5. We have chosen to take all events from M L =1, since we analyze different 
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regions in the field. If the seismometer stations have not changed in this region, all 
events can be taken into account tor a specific region and compared trom year to 
year. 


Table 2-4. The event rate, including Standard deviation, in the regions Central, SW and Other as a 
function of the number of days since the start of seismicity on December 5 ,h 1991. 


Time (days) 

Central 

Event rate 

SW 

Other 

7000 - 8080 

0.058/day ± 

0.017/day ± 

0.098/day ± 


0.006/day 

0.004/day 

0.01/day 

8080 - 8539 

0.024/day ± 

0.039/day ± 

0.091/day ± 


0.007/day 

0.09/day 

0.015/day 


Table 2-3 collects the number of seismic events (magnitude M L > 1) for each 1000 
day period; the corresponding event rate for the 1000 days just before as well as 
539 day elapsed since 17 January 2014 are depicted in Table 2-4. The number of 
events and their event rate is assessed in three regions of the field: “Central” 
(central Loppersum area), “Southwest (SW)” (line Eemskanaal to the area north of 
Hoogezand) and “Other” (the remaining part of the field), see TNO (2014b). 

These results indicate that: 

1) The event rate in the “Central” area has diminished since January 17 th 2014. 

2) The event rate in the “Southwest” area has gone up by a factor of more than two 
after January 17 th 2014 

3) The event rate in the “Other” area has stayed more or less comparable before 
and after January 17 th 2014 

There is no doubt that the Central area experienced a significant drop in seismic 
activity while the Southwest area experienced a significant increase in seismic 
activity. The increase in seismic activity since day 4000 (around 2003) in each of 
the areas is noteworthy, and this is in line with the increase models used in the 
previous report (TNO2014b). 


2.6.3 Seasonality 

In this section the seismic event rate response of the Groningen field is analyzed for 
correlations with seasonal swings in production. 

One way to investigate this is to look at the correlation function between the change 
in production ( dP) on a monthly basis and the number of seismic events ( n ) on a 
monthly basis, Corr (dP, n). Perhaps more tellingly, we may look at the correlation 
between dP and dn - the change in seismic events on a monthly basis, Corr 
(dP,dn). This correlation function is related to the former. The autocorrelation 
functions of dP and dn provide additional information. 

If we compute the above functions for each year stading in 2003 with a maximum 
time lag of 36 months the following results ensue: 

1) The autocorrelation function for dP shows a seasonal trend (Figure 2-10). 
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2) Since 2005 most years show a positive correlation between production 
changes and number of monthly events for 2-8 months (Figure 2-11). 

3) In all years the correlation between production changes and subsequent 
changes in seismic events is maximum and positive at a lag of some two 
months (Figure 2-12). 


The correlations have been evaluated from 2003 since the seismicity over the field 
has been more or less constant up to 2003 and increases after 2003 (Figure 2-9). 
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Figure 2-9. Number of events occurring within the contour of the Groningen gas field as a function 
of time and Magnitude (M). 
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Figure 2-10. Autocorrelation of the production on a monthly basis. 
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Figure 2-11. The correlation between the production on a monthly basis and the number of seismic 
events. 
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Figure 2-12. The correlation between the production on a monthly basis and the change in seismic 
events on a monthly basis. 
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lf we compute the (auto) correlation functions over the whole period 2003 to 2014 
(maximum lag being 36 months again) we note the following: 

1) The autocorrelation function tor dP shows seasonal effects, as to be 
expected (Figure 2-13). 

2) A yearly pattern is obvious in the correlation between production changes 
and number of seismic events: After 5-7 months, 17-19 months, and 29-31 
months the correlation is maximum positive (Figure 2-14). 

3) The two-months effect in 3) is still present, but it is no longer predominant. 
Stacking of all data apparently washes this effect away (Figure 2-14). 



Figure 2-13. Autocorrelation of the production on a monthly basis. 




Figure 2-14a.) The correlation between the production on a monthly basis and the number of 

seismic events and b.) The correlation between the production on a monthly basis and 
the change in seismic events on a monthly basis. 


All together we infer seasonal effects in the seismicity. Since the production 
changes follow more or less identical patterns each year it is not possible to 
attribute the values of the correlation functions at year Y exclusively to the 
production changes in that year. It is well worth remembering that correlation does 
not prove a causal relation. However, what we see is seasonality, whatever the 
precise mechanical processes in the subsurface. 
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2.7 Bayesian point change model 

2.7.1 Introduction 

Change point models are used to detect changes in occurrence rates of events. 
Gupta and Baker (2015) have developed a Bayesian Point Change model which 
quantifies the changes in seismicity rates for Oklahoma US. The unknown 
parameters in this model are the date of change, the event rate before the change 
and the event rate after the change. In Oklahoma a marked increase in seismicity 
was observed after 2008 (Gupta and Baker 2015) and it was confirmed that the 
change point occurs around 2008-2010. Furthermore the post change date 
seismicity rate is 300 times the pre change date seismicity rate, indicating a 
significant increase in the number of seismic events. 


2.7.2 Results of Bayesian Change Point Model for M L > 1.5 and 1991-2015 

In this section the Bayesian Point Change Model (Appendix B) is applied to the 
observed seismic events of the Groningen field. From the seismicity database of the 
KNMI (www.knmi.nl), only the events with magnitudes largerthan M L =1.5 that occur 
within the contours of the Groningen gas field were selected. The magnitude M L of 
1.5 (magnitude of completeness) is chosen as it represents the events which can 
be recorded over the entire field since January 1996. For the analysis a point in the 
center of the field (latitude=53.23 and longitude =6.716) and a radius of 50 km is 
defined, such that all induced events related to the Groningen gas field are selected 
for the analysis. 

Figure 2-15 shows the results of the Bayesian Point Change analysis; a change 
point is observed for January 12 th 2003. In Figure 2-16 the pre change date and 
post change date event rates are shown. The pre change event rate is around 0.01 
events per day, corresponding to 3-5 events per year. The post change event rate 
approximates 0.05 events per day, corresponding to 15-20 events per year. This 
corresponds quite well to the observed seismic events (Figure 2-9). 

After 2009 no change point is observed for all events with magnitudes larger than 
M l =1 .5. The analysis cannot provide a change point for the period after January 
2014 when production was changed over the field due to the relatively small 
number of events after January 2014. 
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Figure 2-15. The probability of change in time over the period of 1991 up to now. 



Figure 2-16a.) The pre change date event rate (in events/day) and b.) the post change date event 
rate (in events/day) 


2.7.3 Results of Bayesian Change Point Model for various locations over Groningen field 

For this analysis an array of points located at approximately 5 km from each other, 
distributed regularly over the Groningen field is investigated. For each location a 
radius of 10 km around the investigated point is taken into consideration. In this way 
the Groningen field is divided in 50 overlapping regions. The probability of change 
in event rate is calculated for each local region and for all events with magnitudes 
largerthan M L =1.5 over the period from 1991 to now. The result is presented in 
Figure 2-17. 

The Bayesian change point model has detected when event rates have changed 
over the entire Groningen field. The earliest change of event rate happens in the 
central part of the field (January 2003). In time the event rate changes spread 
towards the edges of the field. This corresponds to earlier observations of the 
spread of events in time (e.g. NAM 2013). At the south and north edges no change 
of event rate could be detected due to the few recorded events in the 10 km radius 
from the investigated points. 













Figure 2-17. Time of event rate changes evaluated at 50 local points in the Groningen field. 
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3 Compaction field up to 2013 from Inversion 


3.1 Introduction 

In TNO (2013) and TNO (2014b) possible inconsistencies were identified in the 
geological model of the field, mainly by the mismatch between the modeled 
subsidence and the measured subsidence. In these reports we have used the so- 
called forward method, illustrated by Figure 3-1. In the forward model gas 
production is used to model the reduction of pressures in the field. The reduction of 
pressure gives compaction in the field, using a compaction model. Using a transfer 
function compaction can be translated to subsidence at the surface (e.g. Van Opstal 
1974). This forward procedure is sensitive to the quality of the geological model and 
the reservoir dynamical model. As is described in TNO (2013, 2014b), mismatches 
between modeled and measured subsidence were identified leading to possible 
inconsistencies in the porosity and aquifer activity. 



Figure 3-1. Schematics of the forward method. 

The opposite of the forward method is given by the inverse method. In this method 
the measured subsidence is used to compute compaction. The inverse method is 
sensitive to the quality of the subsidence measurements but not sensitive to the 
quality of the geological and the reservoir dynamical model. The identified problems 
in the geological and reservoir dynamical model have led to the implementation of 
the inverse method, described in section 3.2, to provide an alternative compaction 
field for Groningen. 


3.2 The inverse model & doublé differences 

As the inverse method is sensitive to the quality of the subsidence measurements 
the doublé differences measured between optical levelling points have been used. 
Therefore problems with reference points have been avoided. The inverse method 
is described in detail in appendix A in the form of a conference paper submitted on 
May 1 st 2015 to NISOLS (Ninth International Symposium on Land Subsidence). 
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In this study the compaction field from TNO (2014b) has been used as input to the 
inversion. This ensures a new corrected compaction field with similar spatial 
resolution as the previous compaction field. In the inversion procedure correction 
factors to the previous compaction field have been sought. These correction factors 
are applied to the prior compaction field to arrivé at the estimated compaction field. 
Figure 3-2 shows the prior compaction field and the estimated compaction field in 
1993 and 2013. 
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Figure 3-2 Prior compaction fields (top row) and estimated compaction fields (bottom row) in 1993 
(left) and in 2013 (right) 


Figure 3-3 shows the compaction field in 2013 derived from the inversion of 
subsidence data together with the faults in the geological model (NAM, 2013). 
Compaction differs from the prior model, described in TNO (2013, 2014a, 2014b). 
There are four areas (Figure 3-2, bottom right) showing higher compaction (>30 cm) 
of which one is located in the northeast, one in the middle, one in the east (close to 
Appingedam) and one in the south (north of Hoogezand) of the field. Compared to 
TNO (2013, 2014a, b) the area of maximum compaction has shifted to the west. 
Areas characterized by high compaction also seem constrained by faults systems, 
which leads to enhanced differential compaction across faults (Figure 3-3). 

Figure 3-4 shows the same compaction results from inversion in 2013, but now 
including the locations of observed seismic events. The seismic events are 
concentrated in a band from Northwest to Southeast. Contrary to earlier results, the 
areas with high seismic event densities do not correlate with areas of high 
compaction; they correlate with a concentration of faults. This indicates that the 
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faults play a major role in the distribution of seismicity. There seems to be no 
correlation with the offset of the faults (Figure 3-5). Finally Figure 3-6 shows the 
difference in compaction between the model of TNO (2014a) and the resulting 
compaction field derived through inversion. The areas with the largest differences 
correspond to the areas where subsidence was poorly matched (Figure 3-7). These 
areas correspond to areas where TNO has discussed NAM’s porosity estimations 
or areas where active aquifers are assumed in the subsurface model, as has been 
described in TNO (2013, 2014b). 



1:250000 


Figure 3-3. Compaction (m) in 2013 obtained through inversion of subsidence measurements 

(section 3.1). The red line gives the contour of the Groningen field and the black lines 
are the faults that are present in the geological model in Petrel (NAM, 2013). 
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Figure 3-4. Compaction (m) in 2013 obtained through inversion of subsidence measurements 

(section 3.1). The red line gives the contour of the Groningen field and the black lines 
are the faults that are present in the geological model in Petrel (NAM, 2013). Also 
shown is the seismicity in the field, the size of the symbols indicates the magnitudes of 
the events. 
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Figure 3-5. Compaction (m) in 2013 obtained through inversion of subsidence measurements 

(section 3.1). The red line gives the contour of the Groningen field and the faults in the 
geological Petrel model (NAM, 2013) are indicated with their offset. 
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Figure 3-6. The difference between the compaction field of TNO (2013, 2014a,b) and the 
compaction field resulting from inverslon (m). 
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Figure 3-7. Figure 5.13 from TNO (2013). Compaction in 2012 calculated with the RTiCM model 
using the subsurface model. 
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4 Main Findings 


In 2014 the total gas production of the Groningen gas field was 42,41 bcm, which is 
less than the imposed production cap of 42,5 bcm. The Loppersum production 
clusters (LRM, PAU, POS, OVS and ZND) produced in total 2,57 bcm in 2014, 
which is below the production cap of 3 bcm. Production varied over the year with 
the majority of gas being produced during the winter months. 

In support of their advice, State Supervision of Mines has requested the following 
additional technical evaluations from TNO-AGE: 

• An update on the seismicity of the Groningen field 

• compaction field based on inversion of subsidence data 

Update on the seismicity of the Groningen field 

The distribution of higher magnitude (M L >2) events occurring since September 2014 
can be explained by 

The distances to producing clusters vs non-producing clusters which can 
explain the events close to Appingedam 

The increase of production at the Ten Post cluster (POS) in December 
2014 

The first explanation would indicate that the effect of reducing production of the 
Loppersum clusters has been partially overruled by production of other clusters 
close to Appingedam. The second explanation would indicate that sudden 
increases in production could lead to a changing pattern of events in time and 
space. The latter statement cannot yet be proven with statistical significance and 
should therefore be further evaluated. An analysis of production and seismic events 
occurring over time could possibly provide further statistical significance. 

The observed density of seismic events from April 2014 to April 2015 is different 
from densities observed during earlier years. Largest seismic event densities are 
concentrated in the Southwest while the center of the field is marked by lower 
densities. In previous years the density of seismic events in the center were highest 
in the Groningen field (see also TNO 2014b). This indicates that the reduction of 
production in the central area has a marked influence on the number of events in 
the same area. Additionally, there is a striking match of the event density in 
2014/2015 to two known fault systems in the field. These fault systems correspond 
to areas in the field where differential compaction, known to be an indicator for the 
occurrence of seismic events, exists (Figure ii). 

Statistical analysis on the number of seismic events indicates that the number of 
events per day in the center of the field has halved since January 2014. The 
southwestern area, however, shows an increase in the number of events per day. 
Similar to the production, the seismic events of the Groningen fields exhibits clear 
seasonality with a lag of some two months between production changes and a 
change in seismic events. 

A Bayesian change point model that has been successfully applied in Oklahoma, 
U.S.A., has also been applied to the Groningen gas field. A change point in the 
center of the field is found in January 2003. Event rates after 2003 have quadrupled 
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compared to the years prior to 2003. This would indicate that the fault system in the 
center of the Groningen field has reached criticality in the beginning of 2003, i.e. 
small changes in stress will lead to seismic events. Over the whole of the field 
change points are identified which vary in time (trom 2003 to 2010); the earliest 
times in the center of the field and later times at the edges of the field. This 
corresponds to the observation that events have started to occur in the center of the 
field and have spread in time over the field. If the change point indicates when a 
fault system becomes critical then this also means that different fault systems have 
become critical at different times. 

After 2009 no change point in the center of the field has been found tor seismic 
events with magnitudes larger than M L =1.5. The number of events (ML>1.5) since 
2014 is probably not enough to show a change point in event rates. Thus this data 
cannot be used to prove statistically significant changes in event rates since the 
production reduction of January 2014. 


Alternative compaction field 

Inversion of subsidence data has provided a correction to the compaction field 
presented in TNO (2013, 2014a,b). The correction is predominantly applied in 
regions where previously erroneous porosity estimations or aquifer activity were 
suspected. The area of maximum compaction has shifted to the west and does not 
correspond to the area of maximum event density in the center of the field. This 
indicates that in this regard the presence of faults is more important for seismicity 
than the compaction itself. Also differential compaction, known to be an indicator for 
the occurrence of seismicity, is visible over faults. 

This leads to the conclusion that the existing seismological model which NAM has 
used in the production plan (NAM, 2013), based on an empirical relation between 
total compaction and the occurrence of events, needs to be updated. As indicated 
in TNO (2013, 2014a,b) the faults in the reservoir play an important role in the 
occurrence of events within the field and therefore they have to be taken into 
account in any future seismological model. 
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A Inversion of double-difference measurements from 
optical levelling for the Groningen gas field 


This appendix has been submitted as conference paper on May Ist to NI SOLS 
(Ninth International Symposium on Land Subsidence). This conference paper will 
be reviewed by external experts and after revision and resubmitting it will be 
presented on NI SOLS, Japan, November 2015. 

Peter A. Fokker and Karin Van Thienen-Visser 

TNO, Utrecht, The Netherlands 

Correspondence to: Peter A. Fokker (peter.fokker@tno.nl) 

Abstract 

Hydrocarbon extraction lead to compaction of the gas reservoir which is visible as 
subsidence on the surface. Subsidence measurements can therefore be used to 
better estimate reservoir parameters. Total subsidence is derived from the result of 
the measurement of height differences between optical benchmarks. The procedure 
from optical height difference measurements to absolute subsidence is an 
inversion, and the result is often used as an input for consequent inversions on the 
reservoir. We have used the difference measurements directly to invert for 
compaction of the Groningen gas reservoir in the Netherlands. We have used a 
linear inversion exercise to update an already existing reservoir compaction model 
of the field. This procedure yielded areas of increased and decreased levels of 
compaction compared to the existing compaction model in agreement with 
observed discrepancies in porosity and aquifer activity. 

Introduction 

The Groningen gas field is a giant onshore field that has caused substantial 
subsidence since the start of its production in 1963. This subsidence has 
periodically been established by measuring the difference in height of stable 
benchmarks, using optical levelling. Pressures in the field have been closely 
monitored for reservoir management. History matching of the reservoir model on 
the observed pressures has resulted in a reasonably accurate pressure distribution 
development over the field. 

There are a number of parameters in the relationship between the reservoir 
pressure and the subsidence which are more or less uncertain. The first one is the 
compaction coëfficiënt, being dependent on the rock type and the porosity. There is 
also some uncertainty in the pressure estimates in some regions of the field, 
particularly in the connected aquifers, where pressure measurements are not 
available. 

In the present paper we use the raw leveling difference measurements in 
conjunction with the prior knowledge about the Groningen gas reservoir in order to 
constrain the uncertainties. We employ an inverse algorithm to this end, but, instead 
of using interpreted heights, we use the originally measured height differences. In 
an earlier paper we reported the benefits this approach [Fokker and Van Thienen- 
Visser, 2015]. 
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Available data 

The Groningen gas field has been in production since 1963. It is located onshore in 
the Northeast of the Netherlands. Extensive geological, geophysical and reservoir 
engineering data have been used to history-match the reservoir characteristics like 
geometry, porosity and permeability. We had access to the simulated pressure field 
at yearly dates from 1/1/1964 to 1/1/2017. The delta pressures were multiplied by 
the height and the estimated compaction coëfficiënt for each grid cell, based on 
lithology, pressure depletion and porosity. For each x-y location these numbers 
were accumulated over the reservoir layers in order to yield a prior estimate for the 
compaction grid at 9070 x-y locations for 54 times [Van Thienen-Visser et al., 2015]. 
We remapped the provided compaction values to locations on a regular 400x400 
m 2 grid for later manipulation. A map of the input compaction grid and the outline of 
the Groningen gas field in 2012 is provided in Fig. 1. 

In the present study we focused on the use of data acquired through optical 
levelling. Usually, investigators use differences of the interpolated height maps to 
estimate surface movement. The procedure to obtain these differences includes the 
coupling to a reference benchmark or a set of reference benchmarks which are 
supposed to be stable, by integrating along the path of measurements to the stable 
benchmark. This procedure is sensitive to errors in the network and it accumulates 
the inaccuracy of all the measurements in the connecting path. The latter drawback 
can be addressed by providing the full covariance matrix of the resulting height 
estimates; this is, however, rarely done. Also reference benchmarks which, in 
hindsight, are not stable give rise to further inaccuracies. We have therefore chosen 
to use height difference measurements directly. The procedure to obtain double- 
difference estimates has been outlined in an earlier paper [Fokker and Van 
Thienen-Visser, 2015]; it involves the determination of height differences between 
corresponding benchmark pairs in subsequent measurement campaigns, which 
have not necessarily been achieved in the same order. 

Optical levelling campaigns have been performed many times in Groningen with 
different coverage. We had access to a total of 92 campaigns, dating from 1938 to 
2012. Within a total of 7995 benchmarks, more than 26,000 height differences had 
been measured. In this set, 1572 benchmarks had been identified as stable ones in 
the resulting optical levelling database. We have constructed differences between 
stable benchmarks only, using the measurement paths along the unstable ones, 
and used these to construct the doublé differences. Further, we discarded 
benchmarks west of the line with x = 230,000 m and south of the line with y = 
575,000 m in the local coordinate system (RD) to exclude the influence of other 
sources of compaction in those areas (e.g. the depletion of the Annerveen gas field 
south of Groningen). Still, a total of 10860 doublé differences could be constructed 
between 987 benchmarks. The locations of these benchmarks are shown in Figure 
1 . 

Forward model 

Gas production causes reservoir compaction, which, in turn, results in surface 
movement. Compaction in the reservoir may also change certain reservoir 
parameters. For the current study, a one-way coupling suffices - the change in 
porosity due to compaction only affects the reservoir pressure negligibly. We 
employed a linear-elastic model for the subsurface response, with the compacting 
blocks in the reservoir as source terms [Fokker and Orlic, 2006]. Using an influence 
function approach, the subsidence at any surface point then is a superposition of 
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the contributions of all compacting reservoir blocks. For the elastic profile in the 
subsurface we used a homogeneous elastic modulus down to a rigid basement at a 
depth of 5000 m. The reservoir is located at a depth of 3000 m. The connection to 
the doublé differences measured with the optical leveling can be made by making 
the appropriate time differences combined with space differences. 

The goal of the present study was to employ an inverse algorithm on the interpreted 
doublé differences to improve the history match of the reservoir model and the 
predictive capability of the model in terms of reservoir pressures and subsidence. 
We considered the compaction of the reservoir as the uncertain parameter - the 
reservoir pressures and the porosities underlying it would involve too large 
computational efforts for this assessment. To map the uncertainty of the reservoir 
compaction we employed a field of multiplication factors at a spacing of 3200 m in 
space and 4 years in time. Values at the actual grid and intermediate times were 
obtained by interpolation. The prior multiplication values were defined as a constant 
value of unity over the field. A Standard deviation of 0.3 was assumed. The 
mathematics of development 

Inverse model 

For the inverse model we define the vector m as the collection of adjustable model 
parameters, the vector d as the collection of double-difference data, and the matrix 
G, working on the model parameters, as the forward model. The inverse problem is 
then formulated as the task of estimating the vector m for which Gm approaches 
the data vector d best. With additional information present in the form of a prior 
model (m 0 ) and covariance matrices of the measurements (C d ) and of the prior 
model (C m ), the conventional least-squares solution is obtained by maximizing the 
objective function ] given by Tarantola [2005] (or by minimizing —log[/]): 

/ = exp [-i(m - m 0 ) 7 'C- 1 (m - m 0 ) - ±(d - Gm) T C^(d - Gm)] 

For the linear problem at hand, the estimate and its covariance are given by 
m = m 0 + C )n G r (GC )n G 7 ’ + C d ) -1 (d - Gm 0 ) 

= m 0 + (G 7 ’C^ 1 G + Cm 1 ) _1 G r C^ 1 (d - Gm 0 ) 
c rn = c rn - C m G T (GC m G T + C d ) _1 GC m 
= (G%- 1 G + C- 1 )- 1 

in which the first or second line of both expressions can be chosen according to the 
number of data points and model parameters [Tarantola, 2005]. A smoothness 
constraint was added by extending the data vector with a number of elements equal 
to the number of multipliers in the model parameters, and by assigning the 
Laplacian working on m as the forward model for those elements. Furthermore, an 
independent constant vertical velocity for every benchmark was used as an 
additional unknown parameter to allow for movement not caused by the depletion of 
the gas field. 

Results 

The inversion exercise yielded an update of the fields of multiplication values and 
values for the autonomous movement of the benchmarks. With the original unit 
values and with the expected values of the multiplication factors, the forward model 
was rerun. Figure 2 shows the prior and posterior calculated doublé differences 
against the measured values. Although the scatter is still large, there is a clear 
improvement. The quality of the fit, indicated by x 2 = ^(Gm _ d ) 2 /erJ, improved 
from 8.8 to 5.9 - the first and second number being calculated with the prior and the 
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estimated model parameters, respectively. The remaining value around 6, much 
largerthan an optimal value around 1, is presumably related to a remaining 
instability in the selected benchmarks, however it could also mean that the Standard 
deviation of the height difference is too optimistic. The average of the background 
movement of the benchmark is zero; the Standard deviation is 0.5 mm/year. 

There is a clear effect on the compaction fields. Examples of prior and updated 
compaction fields are given in Fig. 3. They show that around some areas the 
compaction levels must be adjusted to explain the measurements. These areas 
consistently return, independent of variations of the amount of smoothing or the 
precise form of the influence function in the forward model. More compaction than 
assumed in the prior model seems to have taken place around Ten Boer [(x RD ; y RD ) 
= (243,000; 588,000)]; less around Delfzijl [(x RD ; y RD ) = (255,000; 592,000) and less 
around Uithuizen [(x RD ; y RD ) = (245,000; 605,000)]. The improvement of the doublé 
difference estimates and the effect on the subsidence estimates benchmarks is 
represented in Figs. 4 and 5. 

Discussion 

The correlation between measured and predicted doublé differences is much better 
for the estimated values of the multiplication factors than for the prior values. Still, 
the scatter remains considerable and there are many points with estimated value 
around zero that show comparatively large measured doublé differences. In view of 
this, it is remarkable that the inversion results in a consistent increase of 
compaction around Ten Boer and consistent decreases around Delfzijl and 
Uithuizen. This result was even apparent when no background movement was 
taken into account and the resulting correlation between measured and predicted 
doublé difference values was even worse. We assume that instabilities of individual 
benchmarks will cause deviations of doublé differences connected to them which 
are compensated with deviations with opposite sign for doublé differences stading 
from them. 

Independent support for the updated compaction field has been found in a separate 
study [Van Thienen-Visser and Breunese, 2015]. In that study, a different forward 
compaction model was employed and the predicted surface subsidence was 
compared to differences of interpreted heights at stable benchmarks and PS-InSAR 
measurement of the surface movement velocity. The areas that we found here were 
also identified in that study, and an additional effort was already recommended 
there to improve the subsurface model in those areas as it pointed towards 
inaccuracies of the porosity model and the assumed aquifer activity. 

Conclusions 

The present study proves the possibility of using doublé differences of optical 
levelling between stable benchmarks for the determination of reservoir parameters 
by its application on the Groningen gas field. The inverse study that we performed 
yielded a consistent update of the compaction of Groningen gas field during the 
lifetime of the field. The area around Ten Boer is compacting more than in the prior 
compaction model; the areas around Delfzijl and Uithuizen less. This is consistent 
with independent results obtained from comparing predicted subsidence with 
temporal differences of interpreted benchmark elevations. A renewed effort of 
reservoir modelling is required to improve the understanding of the reservoir in 
these areas. 



TNO report | TNO 2015 R10755 


Appendix A | 5/7 


References 

Fokker, P. A. and Van Thienen-Visser, K.:. On the use of doublé differences in 
inversion of surface movement measurements. Paper ARMA 15-096, presented at 
the 49 th US Rock Mechanics / Geomechanics Symposium. San Francisco, CA, 
USA, 28 June - 1 July 2015. 

Tarantola: Inverse Problem Theory and Methods for Model Parameter Estimation. 
SIAM, Paris, France, 2005 

Van Thienen-Visser, K and Breunese, J.N.: Induced seismicity of the Groningen 
gas field: history and recent developments. The Leading Edge, special issue 
Injection Induced Seismicity, 2015, in press. 

Van Thienen-Visser, K., Pruiksma, J and Breunese, J.N.: Compaction and 
Subsidence of the Groningen gas field in the Netherlands. Submitted to NISOLS, 
2015 



TNO report | TNO 2015 R10755 


Appendix A | 6/7 


x 10 5 Prior compaction in 2012 



Figure 1 Prior estimate of the compaction field of the Groningen gas field in 2012 
(color-coded), outline of the gas-bearing layers (solid line) and surface locations of 
the benchmarks used in the study (filled dots). 
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Figure 2 Predicted versus measured doublé differences, predicted with prior 
compaction field (left) and with estimated compaction field allowing point noise 
(right). 
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Figure 3 Prior compaction fields (top row) and estimated compaction fields (bottom 
row) in 1993 (left) and in 2013 (right) 
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Figure 4 Improvement of the fit of doublé differences measured for two out of the 
92 campaigns - towards 1993 and 2013 (stading times are variable for the different 
points). The color code indicates the ratio between prior and posterior offset: 

(Gm £ - d)/(Gm 0 - d). Absolute values of this number smaller than 1 (yellow or 
green) indicate improvement. 
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B Introduction theory to Bayesian Point Change Model 


In this appendix the Bayesian Point Change Model is explained. The Bayesian 
Point Change model takes a Bayesian approach to the retrospective analysis of a 
Poisson process with a single change point at an unknown time. The rate of 
occurrence at time s, a(s), is equal to a-i if 0 < s < t and a 2 if s > r The analysis is 
based on the observation period [0,7], during which n events occur at times 
t = (ti.h, ■■■ - f n)- The variables r, a ± , a 2 represent the date of change, event 
occurrence rate before the change, and occurrence rate after the change, 
respectively. 

We assume thatT,^,^ are independent a priori, and that the prior densities of a ± 
and a 2 have the conjugate form: 

p(ctj) oc 1 e aj/e i 

The likelihood is 

£(r, a 1 ,a 2 l t ) = J J cci J”” | a 2 e _ “ 2t 

t=t 1 t=t T +1 

= af (T) e-™i a N ( 7’)-N(T) e -(7’-T)a 2 

where Af(t) is the number of events that occurred in the interval [0, t], Thus, 
posterior density can be calculated as 

p^T.a^a^t) oc £(T,a 1 ,a 2 \t)p(a 1 )p(a 2 ^)p(T^ (3), 

since all parameters are mutually independent. 

The marginal distribution for each of t, a t , and a 2 can thus be obtained by 
integrating the posterior density over the remaining two variables. The posterior 
density of r is thus 


p(r|t) = p(r) f a^ (T:)+kl 1 e ( t+ öi)“ 1 da 1 

ƒ“ a «W-»W+k 2 -i e -(r-T+è)« 2 d ai 

_ 1 r(r 1 (r')) r(r 2 (T)) 

_ T Si(T) ri W S 2 W , 2(:t) 

where ^(r) = Af(r) + /q 
5i(T) = T 4“ ~ 
r 2 (r) = Af (7) — Af (r) + k 2 
S 2 (j') — T — t + y 
Furthermore, 

1 

PW = ^ , 0 < t < 7 

since the prior distribution for the time of change r is assumed to be uniformly 
distributed over the observation period. This means that the change in event rate is 
equally likely to occur at any time during that period. 
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The posterior density of the event rate before the change % is obtained by 
integrating equation 3 over t and a 2 . This does not yield a simple analytic form. The 
function is discontinuous in r. The most convenient form tor numerical integration is 
the sum of integrals of some continuous functions. Therefore, the time range is 
discretized on a daily basis, and summed to approximate the marginal posterior 
distribution. 

T 

pOik) * [i 

T = 0 


The posterior distribution of event rate after the change point, can be calculated in a 
similar way: 


p(« 2 |t) « ^[^ a 2 z(T) le “ 2 - S2(T) r(r 1 (r))S 1 (tY 1 (t) 

T —0 


Bayes factor 

The test for a change point compares a model with a change point to a model with a 
constant event rate. The change point model is applied to the observed data 
assuming that there is a change point. It calculates the probability of change on any 
given date. To actually check whether the data support the presence of a change 
rate or favours the model with the constant event rate, a Bayes factor is used (see 
also TNO 2014b). 

The Bayes factor /? is defined as the ratio of the likelihood function for a constant 
rate model H 0 to that of a change model H t . The constant rate model has only one 
unknown parameter, which is a constant rate of occurrence (i.e., constant seismic 
event rate). For conjugate priory of constant rate the same gamma distribution is 
used: 

p(a 0 ) oc a ; k ° e~ a ° /e o 

leading to likelihood function: 

f co 

£(tf 0 |t)=l T(a 0 |t) p(a 0 )da 0 
Jo 

and similarly for the change model to: 

£(tfi|t) = f g T ƒ” J 0 c °£(r,a 1> a 2 |t) p(a 1 )p(a 2 )p(r) da ± da 2 dz 

lf the value of parameters for gamma conjugate priors are kj = 0.5 and dj -> °° for 
j = 0,1,2 then it is shown by Raftery and Akman (1986) that the equation for Bayes 
factor can be simplified to: 

Pit) = 4V^ r n T(n + 1/2) [X=oF(r 1 (T))5 1 (T)- r iWr(r 2 (T))5 2 (T)-^w] 1 

When the Bayes factor is small enough (less than 1) it means that the change point 
model is supported by the data. In this study, a change point model is favoured if 
the Bayes factor is smaller than 0.001. Every time a change point is foundthe data 
strongly support a change point model, lf, on the other hand, a change point is not 
found (Bayes factor > 0.001); than this means that the constant model is preferred. 
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In the case of a small number of events, the constant model is automatically 
preferred above the change point model. 
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